from kuai.sim import setup_model, SimulationJob
from kuai.mol.io import read_molfile
from kuai.sim.dff2kuai import read_ppffile
from kuai.sim.efunc import make_efunc
from kuai.unit import default_units
from numpy import dot, reshape, float64, matrix, diag
from math import sqrt
from scipy.linalg import eigvals


mol = read_molfile('C2H5OH.msd')

index, parameters = read_ppffile('C2H5OH.ppf')
model = setup_model(mol, index, parameters)
masses = []
for i in model.masses:
    masses += [1/sqrt(i)] * 3
    
masses = diag(masses)

funcs = make_efunc(model, index)

# for func in funcs:
#     print func.name, len(func.lists)

fe, _, _ = default_units.format(1, 'kcal/mol')
fg, _, _ = default_units.format(1, 'kcal/mol A')
fh, unitU, unitO = default_units.format(1, 'kcal/mol A A')
fHz, _, _  = default_units.format(1, 'Hz')

print unitU
print unitO

job = SimulationJob()
e0 = job.getE(model, funcs)
e1, g1 = job.getEG(model, funcs)
e2, g2, h2 = job.getEGH(model, funcs)

h2 = matrix(reshape(h2, (len(model.coords), len(model.coords))), float64)

"""
print (e0 - e1) * fe
print (e2 - e1) * fe
print (g2 - g1) * fg
print e2 * fe
print g2 * fg
print h2 * fh
"""

h2 *=  fh*11791.0
v = masses*h2*masses

for i in eigvals(v):
    i = i**0.5
    print abs(i), i
    # print i
    # print i * 5308.4946858738958

"""
for i in range(2, 16):
    job.setCPU(i)
    e2, g2 = job.getEG(model, funcs)
    deltae = e2-e1
    deltae = sqrt(dot(deltae, deltae))
    deltag = g2-g1
    deltag = sqrt(dot(deltag, deltag))
    print i, deltae*fe, deltae/sum(e1) * 100, deltag*fg
"""
